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ABSTRACT 
Finite difference eS to the one-dimensional heat 
- 
equation U, = U = are used to introduce 6Xplicit and implicit dif-- 


ference equations, The convergence and stability of these equations 
is discussed and these concepts are used to show the restrictions 
imposed on explicit equations. The Implicit Alternating Direction 
(IAD) method is introduced as a means for solution of the two-di- 
mensional heat equation, Although the IAD method in its basic form 
is applicable only to parabolic problems, it is possible by slight 
modification to apply the method to elliptic problems, Two examples 
are used to illustrate the use of the IAD method for solution of 
parabolic and elliptic equations for a rectangular region, These 
examples include a work requirement comparison with other difference 
methods, A third example is given to show the applicability of the 
IAD method to non-rectangular regions, in this case a parabolic 
problem over a circular region, Results of these examples show that 
the IAD method is a convenient and accurate method when applied to 
both parabolic and elliptic partial differential equations and 


suggest applicability to a wide range of problems, 
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INTRODUCTION 

U.+U. = cu, is the parabolic partial differential equation 
governing the flow of heat in two space dimensions, The constant k 
appearing in this equation is the diffusivity of the material through 
which heat flows, We assume here and subsequently that the time unit 
is chosen such that k = 1, Numerical approximations to the solution 
of this equation may be obtained by solving an associated difference 
equation, Depending upon the time level at which the derivatives 
Ue and Uy are approximated, the difference equation will be one 
of two types: (1) explicit, which is simple to solve since it involves 
only one unknown quantity in terms of several lmown quantities, but 
which may also require an enormous number of calculations due to 
limitations on the size of the time-step; and (2) implicit, which 
does not restrict the size of the time-step, but which involves five 
unknown quantities and is most easily solved by an iteration process, 

A modified implicit difference method, called the Implicit 
Alternating Direction (IAD) method, was introduced by Peaceman and 
Rachford in 1955 and may be used to obtain numerical approximations 
to the solution of the heat flow equation in a direct and convenient 
manner, in addition, the IAD method may be applied, with very slight 
modification in form, to elliptic partial differential equations, 
The JAD method involves two difference equations: the first in 
which three unknown values of U in the x direction are expressed in 
terms of three known values of U in the y direction; and the second 
in which three unknown values of U in the y direction are expressed 
in terms of three lnown values of U in the x direction, Each of these 


equations gives rise to a system of simultaneous equations, However, 


the coefficient matrix of each system is tri-diagonal, and therefore 
the solution is relatively easy to obtain, 

In this paper we will first consider the simple one-dimensional 
heat equation us = U,. We will derive the associated explicit and 
implicit difference equations and then introduce the concepts of 
convergence and stability to show that a restriction on the time- 
step for the explicit equation is necessary, The explicit and im- 
plicit difference equations for the two-dimensional heat flow equa- 
tion will be given and the associated restrictions and solution 
difficulties will be considered, 

We shall then examine the difference equations which define the 
IAD method, showing how these equations are solved and how they are 
applied to the elliptic problem governed by Laplace's equation 


U U_ = 0, We will discuss some aspects of the proof that the 


xx * yy 

IAD method is both convergent and stable, particularly in regard to 
selection of an iteration parameter for use in solving elliptic 
equations, 

Three examples are given, The first example is an application 
of the IAD method to the parabolic partial differential equation in 
two space variables, and the second is an application to the related 
steady state problem, The region for the problems in both of these 
examples is rectangular, The third example is an application of the 
IAD method to the same parabolic equation, written in polar coordi- 
nates, over a circular region, In Examples one and two a compar- 


ison is made of the work requirement for the IAD method and the work 


requirement for other difference methods, Data obtained in Examples 


Ad 


one and three are used in an attempt to obtain bounds for the 


discretization error, 
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I, Comparison of Explicit and Implicit Difference Equations for 
the Problem U, = U,,. 


For simplicity, consider the parabolic partial differential 
equation (PDE) 


(oD ene x,t) = UL (x,t) for a<x<b, t)d0 


. 
with bsandary and initial conditions 

Wane) = et), UCb,t) er h(t), Uix,0) = f(x). 
This is the equation that specifies the conduction of heat in one 
space dimension, 

To attempt the solution of this problem by a finite difference 
method, we establish a network of grid points throughout the region 
ae t>0, In the x direction divide the region into strips of 
width x = (bea)/M, and in the t direction into strips of width 


At = T/N, where M and N are arbitrary positive integers, and T is 


some fixed time for which the solution is desired. We may now use 


indices i and k to denote the grid points 14x and kAt, 
If we assume that U(x,t) possesses a sufficient number of 
partial derivatives, then the relation between U(x,t) and 


U(x +ax, t +4t) is given by the Taylor's series expansion 


U(xrax,t+at)= Ulx,t) + [ax 2 tat 2 JuUMt) + 


2 2 mi 
+, [ax 3 rtd] Uwt)r °° tt EX3 t OtH| U(r) + Ry ) 


where the remainder term is 


n 
R = (1/nt), ax do + ot de ]"u(K+ vax, t +Xat) 
for 0£%£1, 0f€X£1, This leads to the following definition: 
Definition, If, in the relation between U(x,t) and U(x + Ax, 


t + At), we neglect the remainder Rw then we introduce an error 
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(called the truncation error) which is of order ( Ax re) ae denoted 
O( ax + at)”, By this we mean that there exists a positive constant 
C such that 
[R| = C( Ax + At)” as both ax and At go to zero, 
Now if we let Vane be the finite difference approximation to 
U(x,t) at the grid point (i 4x, k At), then by using Taylor's 
expansion we obtain the following difference approximations for the 


derivatives U, and U 
t xX 


V V5 / ot + O0( At) 


t~ (Vs iced - 


oo daeaae eS 


2 2 
~ ov, +My /(ax)® + 0( ax)*, 
Using these approximations, one finite difference approximation to 
the PDE is 
Z 
2 a es Vs el DEAT, oa Viet uo/(42)"- 


If we define Y = ( ax\/at and then solve for v we get the 
41 ,k+1 


equation 
2) Vs nea =O) g ie + H2/P ) vy HOP) Mays 
The boundary and initial conditions are approximated by 


YO +d =g [ (k+l) ot] ; YM k+l = h[ (k+1) ot | -, and ~ 


Ys 0 = r(i Ax), where i= 0 and i = M represent the points x = a 
and X = b respectively, and k = 0 represents t = 0. Equation (2) 
is an explicit difference equation since we can calculate the values 
of er explicitly at time level (k + 1)at if we know the values 


of we at time kA4t, Thus by starting with the boundary and initial 
t 
conditions, equation (2) can be used to calculate the values of v 


at all grid points for any time t 90. 


L5 


If we approximate the derivative U, at time level (K+ eat 


rather than at time level kat, we obtain the difference approximation 


Sane Vs fat ~ Caieg - et i Vi ore. 


Again, letting Ie = (aox)*/ot and rearranging, we have 


OV ie | sey a ee af, 

Equation (3) is an implicit difference equation since it involves 
three unknown quantities at time level (k+l) 4t in terms of one 
known quantity at time kAt, The boundary and initial conditions 
have the same approximations as for the explicit equation, Thus 
we have a more complicated scheme in this case since at any time 
level (k+1) St, the implicit equation (3) is written once for each 
point 1£i€£M-1l, resulting in the following system of M-1 simul- 


taneous equations in the M-l unknowns La ee 


Cit Ye)Viner 7 Ae Va, key =Vnt % g[(K+ Nae] 
-Vp Viycd H(1+ Be) Va, a = Vo Va Kel = Vz,k 
a 4p Vert, ie eM Wey Ve, w+ ~ tp Viel, Kel = Vi,K 

- /p Ve K+ + (I+ Hoenn, K+ = Vu-i,K * Vp h[(k+ iat] 


Now let a, =c, = “l/p and b, = (1+2/y Ne (i = 1,2, . 9 ), 
and call the know terms on the right hand sides of these equations 
dj» d., re quel? Then after dropping the subscripts k+l, the 


system has the following form 


bv + C1V5 = 
aoVy + boV5 + CoV3 = d., 
my, | + byvy 53 5549 sim *'* dy 7 
ay-3u-3 * Py-o%M-2 + Sm-2%M-1 = Sn-2 
ere "ye2 * wea wwe 
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The matrix of coefficients for this system is tri-diagonal, and 
the use of GausSian elimination results in the following algorithm 


for solution, 
Let 3,= b, ) af = 41/3, 
(= bi- (Ail VB ;, » B= (der ag Mey az , (i=2,3)+° +) MHI) 5 


then Vet ~ {hee ) Viz Ur Ce Miri YB ) Cis Maan M=5, ee Al, e 


Thus we have our choice of two finite difference schemes with 
which we can attempt to find an approximate solution to the PDE (1). 
We shall see that one of the conditions required to make the 
approximation a useful one is a restriction on the size of the time- 
step At when using the explicit equation (2). It will be shown in 
Example 1 that the use of the implicit equation (3) requires less 
work than the use of (2), even with the added complexity of solving 
systems of equatiors since there is no restriction on the size of At 
in the implicit method, 


Whether or not the sylution v of the difference equation is 


tek 
a good approximation to the solution U(x,t) of the PDE can be answered 


by considering the questions of convergence and stability. We begin 
with the following definitions: 


= U(iax, kAt) - v is called 


Definition, The difference w 4k 


a ire 
the discretization error at the point (iax, k4t), 


Definition, We say that the solution v of the difference equation 


i,k 
converges to the solution U(x,t) of the PDE at the point (i 4x, kat) 
if lin W. == 
At X90 15k 

(i ax, kat), then we say that the difference procedure ts:convergent, 


0. Ifv, , converges to U(x,t) at every grid point 


When we speak of convergence, we are speaking about the difference 


between the exact solution to the PDE and the exact solution to the 


1? 


difference equation, Since we will generally attempt to solve the 
difference equation using a digital computer, we cannot expect to 
obtain the exact solution, Instead we obtain a numerical solution 
which differs from the exact solution due to round-off error, When 
this difference remains small, we say that the corresponding procedure 
is stable, 

Before giving a formal definition of stability, we examine 
the question of round-off error in more detail, Suppose v(x,t) is 
the exact solution to the difference equation and suppose we replace 
v(x, t.) by v(x, »t,) +e at the grid point (x, t.). This may be 
considered to be the effect of round-off error or possibly a mistake 
in the calculations, We call e the error at the point (x to). lige 
we then continue the calculation using the value v(x, a) + e, with- 
out introducing new errors, we obtain a solution v*(x,t) at the 
point (x,t). We call the difference v*(x,t) - v(x,t) the departure 
of the solution caused by the error e at (x, t.). If errors are 
committed at more than one point, we speak of the cumulative depar- 
ture due to these errors, For linear problems we can use the 
principle of superposition to conclude that the departure due to two 
or more errors is the sum of the departures due to the individual 
errors, We can now give the following definition: 
Definition, A difference procedure is stable if the cumulative 
departure, due to errors in the solution, tends to zero as the 
maximum of these errors tends to zero and does not "grow faster" than 
some positive power of 1/Ax as 4x tends to zero, We say that the 
cumulative departure does not "grow faster'' than some positive power 


of 1/ax to exclude exponential growth in 1/Ax, 
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The general method used to show convergence is to find a bound 
for the discretization error w(x,t), and then show that this bound 
depends on ax and st and tends to zero as Ax,OSt>0, Stability 
is most often proved by a method due to Von Neumann, This method 
involves expanding the error in a Fourier series and showing by means 
of the coefficients of this series that this error does not grow as 
the solution progresses, For the two difference schemes which we 
have introduced, we will obtain the conditions for convergence and 
stability by a method given in Forsythe and Wasow [3] . 

We assume, for convenience, that the associated boundary condi- 
tions for the PDE (1) are g(t) = 0 and h(t) = 0, Without loss of 
generality we can assume that a » 0 and b =TT since we can modify 
this interval by means of a linear transformation. Then the exact 
solution to (1) is 


2 
(4) U(x,t) = 2. ae \ sin rx, 
=| 


where the initial function f(x) has been expanded in the convergent 
Fourier sine series 
od TY 
A - a 
f(x) = Z- ! asin rx, where a. = 2{ f(x) sin rx dx, 
a (=) 


Now consider the following difference equation 
2 
(5) 5 reg Mic St = Fg 47 e425 lead + Vand yen)/ (OX) 
2 
+ (l-a) (Vs ak e"4 + V5 211) /C ax)”, 
where oO” is a parameter, Note that for & = 0, equation (5) 
reduces to the explicit difference equation (2); and for 7 =1, it 
reduces to the implicit equation (3). We can now examine the stabil- 
ity and convergence of equations (2) and (3). 


Theorem, When o < 1/2, equation (5) represents a convergent and 
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stable difference scheme if a : CePA t => 2(1-20), When 
O 1/2, equation (5) represents a convergent and stable difference 
scheme for all values of 2 5 


Proof. Suppose that v(x,t) is a solution to (5) which satisfies the 





conditions: v(x,0) = f(x), 0<x<mT ; v(0,t) = vG@r ,t) = 0, 0<t ST, 
Here, for convenience of notation, we use (x,t) to mean the grid point 


(iax, kat). Assume that f(x) can be written in the Fourier series 


Mis 


al Ge), a Sin rx, where 
r 


ae f(x) sin rx dx and > | | oe 
v= | 


+ 


we now seek a solution of the form 

v(x,t) = 2. a. Y (t) sin rx, 
if Y' (0) = 1, the solution satisfies the initial condition; the 
boundary conditions are automatically satisfied, If we substitute 
W (+t) sin rx into (5), we get . 


With the initial condition Y (0) = 1, this difference equation has 


the solution V7, (t) = gp 


where @D . S(i-oc)sm' = os 
7 Heer Vax 


and where t/At takes only positive integer values. So we have as 


a solution 


Oo} 
t/At 
(6) Vist) = : a P) sin rx, 
Ve 


if this series converges, Since > [ay <<” , a sufficient 
ele 
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condition for convergence (in fact, for uniform convergence) is 


IP | £ 1 for all r, ‘But this means 
lee Yo (\- ao) sin*v Sf 
i+ 42 sin*y dx 

. z 





eal 


The right hand inequality is always satisfied, The left hand 
inequality is satisfied when 2(1 - 20)< e . Note that this is 
always the case when ©21/2, so that [P| & Peforvany p « 

To show convergence it is necessary to show that the series 
(6) tends to that given by (4) as at, ax— 0, 

; wae 2 2 
In the expression for P. we replace sin“r ax by r-(Ox)"/4 
za 


+ OW) to get 
[-v ‘at +At|o vi+ $(i-o )O(ax)* | 


ic - }+at[ovi+4cO0@xPl 
tht 
We now examine the Lim QD . if we consider this as an 
OX At 0 V 


iterated limit, we have two cases, These are: 


ht I—-vtatt ovat 
Case l, Lim Lim PD = La} [eae 
At9O AX20 'V At-70 [+ Gavi Are 


tKt 


4 aes 
- Lim |i-v7at 4 Nt) | =~ 


AtJ0 [tovVviaAt 
2 
= } Lim fit ot + ZVGe) Vise vee 
At70 lt ovat = 
t/a t | 
Case 2. Lim Lim F As it stands, Lin Dp t/a 1# is 
4X70 AtVO r ey ie 


t = 
; : ; {_ QD pac_ i e* in, 
an indeterminate form so we write im ( = Mm 
L At IO Ato 
tLim 
= eC atvo At 


Using L'Hospital's rule we have 


Lini in fe = t im a Ly Q 


pe 
¢ dim ae ee my = t[40-20)00 x)~v a 


“all 


—£140-20) OG@a)*-V? J 
eC = 


tAt re 
eo {in Lin Y = lela -vt 


AXIO Atva © OK Pe 


Thus we see that regardless of the way in whichax andat tend 
Z 
to zero, Lim Gd t/a t = e § so that every term of (6) 
At sxoo' Pr 
tends to the corresponding term of (4), 

Since the convergence of (6) is uniform, it is possible to 
consider the limit, as Ax,At tend to zero, in a termwise fashion, 
Therefore if 2(1 - 2a) ~ , then lim v(x,t) = U(x,t), and the 

AtAX30 
difference equation (5) represents a convergent scheme. 

To show stability it is more convenient to represent the 

values of v(x,t) at each grid point by means of a finite Fourier 


series rather than the infinite series (6), This is based on the 


theory of trigonometric interpolation and uses the orthogonality 


relations 
<4 Cos “Cos Ne 
° CO os a 
> gan ynvax 40" } mvax =O nezm,nm<M 
Vi=O 
mI M-| 
mt St tmaxe Sh uit 
/ cos nrox= sin” nrAx = M/2, 0<n<M, 
Yeo Va O 


where M = 717 /A x is an integer. Using these relations, we can 


M-! 
f(x) = 2 A. sin nx, 
mo, 


write 


where m= 
AL =(2/M ) = f(r Ax) sin nr Ax, 
This relation will banvalae for the points x = rAx,(r = me Ma). 
we use a finite series since round-off error cannot be regarded as 
a smooth function independent of Ax, 


Now, let M = 7TT/ax be the number of grid points on a line for 


constant t and consider an arbitrary line of errors ej(p=l,2,... M=1) 


ae 


at t = 0, We then consider these errors as some initial function 


and represent them by 
M-| 


e = r ; A. sin sr Ax 
S 1 


Ves 
(s =1,2,...,M-1) 


where M-| 


A. = (2/™) S mos Sane pr Ak 


Cee ey ee 
Then by the same argument that resulted in the infinite series (6), 
we show that ae the 
(7) v¥(sax,nAt) = > ao sin sr Ax, 
vied 
is the solution of (5) at the points of the grid for the initial 
values es» or v* represents the departure due to e. for the linet. = 0. 


We now determine a bound for v*, Substituting the expression 


for Ne in (7) we have 


M-| 
* = 
(8) v*(s 4x,nat) = > ee 
P= 
where <f 
M thet 


Bs ,n - (2/M) pa Y sin sr Ax sin prax 
v= | 


(Here, our interchange of summations is valid since these are 


finite series), If 2(1-20)<€ ie , then Z| < 1 and le, nl = 2, 


Let lOere § for r=1,2,...,M-l where § represents the maximum 
of the errors, Then from (8) we have 
v¥(sax,ndt) ££ MS = 2178 /ax 
which means that the cumulative departure due to errors is of order 
5 ( ax) which by our definition means that the process is stable, 
Therefore, we have proved that for 22(1 -~ 20°), equation (5) 


represents a convergent and stable finite difference scheme, It 
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follows then that for o > 1/2, (5) represents a convergent and 
stable scheme for any positive value of i . Forse =10,ae) 
reduces to equation (2) and the condition for convergence and 
stability is P22, or At é (ax)*/2. As previously stated, the 
explicit method is restrictive in the size of the time-step. For 
Oo = 1, (5) reduces to the implicit equation (3) and is unrestricted 
as to the size of the time-step allowed. 

The method used to prove stability in the previous theorem is 
essentially the method of Von Neumann, This method can be extended 
to problems involving two or more space variables, 

Note that for the difference equations examined, the conditions 
for stability and convergence are identical. Lax and Richtmyer [5] 
have developed a theory for a large class of linear partial 
differential equations with constant coefficients which, for slight- 
ly different definitions of convergence and stability, states that 
stability of a difference equation implies convergence of the equa- 
tion, Briefly, this theory assumes that the initial functions of 
the linear differential problem belong to some Banach space 8B, of 
functions of x, Then associated with every initial function in B 
is a solution to the difference problem (this solution being a 
function of x alone when t is fixed) as well as a so-called genuine 
solution to the PDE, Both of these functions belong to B, Under 
these conditions a finite difference procedure is called convergent 
if the solution to the associated difference equation converges to 
a genuine solution for all initial functions in B, This is a more 
restrictive definition than the one we have used since we have 


established convergence for initial functions in some incomplete 
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subset of the space B, A difference procedure is called stable in 
this sense if the solution v(x,t), corresponding to some initial 
function f(x), satisfies 

Ilv(x,t)]] € mt) {| f(x) I] , 0 € t ¥ T, for all f(x) in B, 
where M(t) is independent of sax, This definition is also more 
restrictive than ours since we allow a bound which depends on Ax, 
There is also a difference in the concept of distance between two 
functions, We have used the norm 

ex, Ceo tyr = Sup [vix,t) = UCx,t)1. 4 

QAtxeb 
whereas Lax uses the mean square norm 
K zl V2. 
ff [v x,t) - UG) dx § 

If the definitions of convergence and stability are taken in the 
above sense, then convergence and stability are proven to be equiv= 
alent concepts in the Equivalence Theorem of Lax, It can be shown 
that explicit difference equation (2), with the restrictions on At, 


actually does satisfy the Lax~Richtmyer theory, 


é? 


Il, The explicit and Implicit Difference Equations for the Problem 
in Two Space Variables, 


For the problem U.. +U_ = U,, we use the same difference 
xX yy ie 
approximations for Ue and U, and a similar approximation for U 


(the approximations for ae and os taken at time level kat). We 
change notation slightly, letting the subscripts i and j refer to 
increments in the x and y directions, and the superscript k refer 
to increments in time, The explicit difference equation is 
(9) an = a + 5 Men, + Very is ay + Ve, ye07 4 Ver | 
Again taking the approximation for _ and Uyy at time level 


(k + 1) At results in the implicit difference equation 
K+ 


(10) vie Ven, V+ Ves - (++ Y) Ves ~is ie Vy 


~t 2 4 Gr 


As in the case of the difference equations in one Space vari- 
able, we find that the explicit precedure (9) is both stable and 
convergent for certain values of P and the implicit procedure (10) 
is stable and convergent for all values of e. (Here we have taken 

p = (ax)*/o he Car ae The restriction on ie for the 
explicit equation is 2 4, This result is proven in Douglas 

[6] » an which the difference equations are treated as matrix 
operators so that the result is valid for any number of space dimen- 
sions, 

For the explicit equation, 24, means that At < (ax)*/4, 
so that the allowable time-step is half the size of that allowed 
for the explicit equation in one space variable, This may mean that 
an enormous number of calculations are required to obtain a solution 


for a particular time T, 
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Equation (10) again results in a system of simultaneous 
equations which can be solved by Gaussian elimination. However, 
such a method for solution is impractical as there are now five 
unknowns per equation and the number of required computations is 
more than twice the number required for the tri-diagonal system, 
The solution of this equation is possible by iteration techniques, 
but we shall demonstrate the solution of the PDE in two space vari- 
ables by means of a modified implicit method called the Implicit 


Alternating Direction (IAD) method, 


a 


III, The IAD Method 

We now introduce the IAD method of Peaceman and Rachford [7] : 
If only one of the second partial derivatives for the PDE, Lean + Uy 
= UL. say Ue is replaced by a difference approximation at time 
(k +1)4t, and a is replaced by an approximation at time level 
kAt, then we again have a difference equation with only three unknowns, 
and this equation may be solved by the use of the previously stated 
algorithm, The difference equation obtained is said to be implicit 
in the x direction, We then repeat the process for a second time- 
step of equal size, reversing the level at which the second deriva- 
tives are approximated, to obtain a difference equation implicit in 
the y direction, When writing these equation, we use time levels 
2k + 1 and 2k + 2 to emphasize that the time-steps must be of equal 
Size, The equations are: 


Implicit in the x direction 


any aK uke 2K+ 2Kt ke LK 2K 
(11) Vey | 2 Ve, = “KY Very + Vert, | + Ven yo “AM, Vey 

4% (4x )* a (ay)* 
Implicit in the y direction 

AK+Z in¢ | , 2K! a ZK+ | 2Kr2 LKt2 2KFL 
ge) ous Moy = ey 7 RM t May Menge 2 Mey Mepoy 

At/y, (ae (a4)* 


We choose Ax = Ay, and let i = (ax)*/aot = ( ay)“/ At, Then, 
rearranging, we obtain a form more convenient for calculation: 

2Kt | yoo Wee ZK 
Qi) am i me ROP +) Vey) = Vet = oe =| T R(P-1)v GJ + UL, ti 


Kee LRLE 2K+ | Z2Kt | Z2Ktl 


(12a) ng) tAC P+), — Vey = Vee ys) ti(y- ik, t Vegi, | 
bach of these equations leads to M ~- 1 equation in M - 1 un- 
knowns at each half time-step, At/2, The systems each have a tri- 


diagonal coefficient matrix, 


28 


In addition to the time-dependent problem U, + U_ = U,, the 


yy 0’ 
IAD method may also be used in the solution of steady state 
problems; for example, Laplace's equation, eile = 0, The pro- 
cedure is essentially one of iteration with each stage of iteration 
being regarded as a time-step in a time-dependent problem, A start- 
ing value is used which approximates the boundary conditions and 
this serves as an initial condition, Equations (lla) and (12a) are 
used as alternate stages of iteration with ig serving as an iteration 
parameter, We will show by means of examples how the IAD method is 
used in the solution of both time-dependent and steady state prob~ 
lems, 

We now state the following theorem without proof. 

Theorem, The IAD method, as applied to the partial differential 
equations AU(x,y,t) = U, (xy ,t) and AU(x,y) = 0, having either 
Dirichlet or Von Neumann boundary conditions, is both convergent 
and stable for any region, 

The proof of this theorem can be found in Forsythe and Wasow 

[3] » pp 272 and 411, and in Birkhoff and Varga [3] . We now 
state some of the results of the proof given in [3] for the IAD 
method as applied to the steady state problem, 

Consider the difference approximations for UL. and Dy used in 
equations (11) and (12), If we regard the solution of the difference 
equation as a column vector V(x,y), then the symmetric difference 
expression used to approximate Ube can be regarded as a matrix an 
operating on V, Thus, 

AV (x,y) = =V(x=Ax,y) + 2V(x,y) - V(x+4x,y). 


Similarity, Lor Uy? the analogous matrix is AY, and 
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A'V(x,y) = -V(x,y- oy) + 2V(x,y) - V(x,ytoy). 
For the case of Laplace's equation in a square region with ne 
investor grid points, ah will be an nx n* matrix which may be 
partitioned into ann xn block diagonal matrix with each diagonal 
block being ann xn tri-diagonal matrix, AY will be similar to 
such a matrix, If we use the subscript k to denote the iteration 


number, the IAD method is then defined by the equations 


Gs) Vio ig - Vice eZ (A' Vk. y, r AY V1) 


Fc 
(4) VK = Veen” A Ven a AY VK ) 


where the parameter y= aaa may vary from iteration to iteration, 
In effect, we are solving the elliptic (steady state) problem as the 
asymptotic solution of the time-dependent problem; thus the At 
appearing in the expression for A: In the actual solution process 
for the elliptic problem, however, At will not enter into the method 
for determining optimal values of y - 

We now define the error Ev as the difference Vy - U, where U 
is the solution of the PDH at each grid point and may thus be re- 


carded as a vector, Then, after eliminating V between (13) 


K-1/2 
and (14), the effect of one complete iteration is 


by = P( Px )Ey 


where the matrix operator P is 
P(e) = (1+ SAY) (I- fe AY)(I+ SAM) (I> AY) 
and I is the identity matrix, 
Now, as shown by Forsythe and Wasow, a sufficient condition for 
vonvergence is that the spectral radius of P(Y) be less than one, 


where the spectral radius is the maximum of the absolute values 
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of the eigenvalues of P, This is shown to be true for any positive 
constant value of ~: 

If the region of interest for the PDE is a rectangle, spec- 
ific statements can be made about the values of the iteration para- 
meter Px which lead to the most rapid convergence, As shown in 
Birkhoff and Varga [3] » however, this is not possible in non- 
rectangular regions, Suppose, for convenience, that the boundary of 
the region is a square of side lengthtTl., Then oe and ce commute 
and, by a theorem due to Frobenius, both matrices share the eigen- 
vectors 

sin px sin qy, (p,q = 1,2,...,M-l), 
with eigenvalues 
Ly ie p ax/2, 4 sin’ q dy/2, (p,q = 1,2,... Mel). 
The eigenvalues of P(p ) are then 
—~ 4 51n* par = 204 _ 
ae or zero ek a 
From this expression it can be seen that the eigenvalues of P(Y) 


(15) 


are less than one in absolute value for any positive value of ~- 
Peaceman and Rachford [7 | use the Von Neumann method of stability 
analysis to derive the same expression as an amplification factor 
for error, In Example 2 we will show how this expression is used 
to determine the optimal values of the iteration parameter ) . 

We will show in Examples 1 and 2 how the IAD method can be 
used to obtain an approximate solution to the time-dependent and 
the steady state problems for a rectangular region, In Example 
3 the method will be applied to time-dependent problem over a 


section of a circular region, 
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mxample 1, Solution of a Time-Dependent Problem Using the IAD 
Method, 


Problem: To find the temperature distribution for various 


times in a square plate with edge-length one. 


The temperature distribution in a square plate of edge-length 
one (See Figure 24) is governed by the following parabolic partial 
differential equation: 

U (x.y, t) > Uy (xo¥ ot) = UL, &%,y,t) , 
with boundary conditions 

U(l,y,t) = U(x,l,t) = 0 UL (O,y,t) = U(x,0,t) = 0 
and initial condition 


Olesen, ,0) = all 





o Uy =90 
Figure 24 Temperatures in a Square Plate 
To solve this problem using the IAD method we use the differ- 


ence equations 


2K+I int 2K+| ZK LK aK 
cua) - Met, 3 2 (Pte, - Vet} = VA +2 (f-1)%,, + Vere 
2K12 2K+2 Z 2K +i 2K! ZiK+ ( 
(12a) Veja cf 2CPHIM = er = vin +2(p-iM,, t Meaty | 
2 u 
where e = sed = ey) 


aya 


fquation (lla) is said to be Implicit in the x direction and 
equation (12a) is said to be {mplicit in the y direction, 

To apply these equations, we divide the region shown in Figure 
24 into N increments of length 1/N in both the x and y directions, 
Here N is an arbitrary positive integer. The resulting grid is shown 
in Figure 25. Thus wget refers to the difference approximation to 


the solution U(x,y,t) at the grid point (iAx, jAy) at time (2k+1) at , 


where 4x = say =1/N, 


J 







(0,Na4) HON Nag) (lox, Nay) 






. Pa 
(0, ja4) — PED) SF Coley yoy) 


(¢ax, 0) (NAX, O) 
Figure 25 Grid for the Region Shown in Figure 24 


Hquations (lla) and (12a) can be used for all interior points 


of the grid and, of course, the values of U are known on the boundar- 
ies x = 1 and y=1. However, at the boundaries x = 0 and y = 0, the 
boundary conditions are U,(0,y,t) = 0 and U.(x,0,t) = 0, so we must 
find an approximation involving UL. for Vis at x = 0. and an approx- 
imation involving Ly for Moy at y = 0, Since the procedure is the 
same at both boundaries, the details are given for the boundary x = 0 
only. 

We assume, for a fictitious point, a distance ax outside the 


boundary, that U(x=- ax,y,t) = U(x+ ax,y,t). This is intuitively 
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acceptable since ee and the region under consideration can be 
thought of as one quarter of a square plate having edge-length 
two and being centered at the origin. For this larger plate (with 
boundaries held at temperature zero) there would be complete sym- 
metry and hence no heat flow across the x and y axes, We now write 
the Taylor expansions for U(x- Ax,y,t) and U(x+ax,y,t) in terms of 
U(x,y,t), where (x,y,t) represents the boundary point 
[o,jay,(2k+1) at]. Thus we have 


' = © (ax)* 
U(x- ax,y,t) UG) OxU, (x,y,t) + ai Uf (x»¥»t) 


3 
= ot) U(x t) + cece 


and : 
U(x+ ax,y,t) = U(x,y,t) + AxU_(x,y,t) + GX U_(x,y,t) 
x Zz! xx 


(ax) 
ay SH Uber st) 8 oe 8 @ 


Now, adding these equations, and neglecting terms involving U 
and all higher derivatives, we get 


U(x= ax,y,t) + U(x+ Ox,y,t) = 2U(x,y,t) + (ax)°U (x.y st) + 0( 4 a 


But U(x- ax,y,t) = U(x+ax,y,t), so we can write 


(16) UL (xy,t) = Zs [Ulxt a xayt) - U(x.y.t)| + 0( ax)*, 
(17) UL Gxyst) = rer A ay,t) - U(xy.t)| + 0Cay)*. 


For notation in the computer program IADMT, we use TST(I,J) to 
denote the values of our approximation to U(x,y,t) at the end of a 
first half time-step (using equation lla) and TEND(I,J) to denote 
the approximation at the end of a second half time-step (using equa- 
tion 12a). Thus TEND(I,J) represents the numerical solution at the 
end of a complete time-step, In the equations that follow, we 


abbreviate TST(I,J) by T*(I,J) and TEND(I,J) by T(I,J). 
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All computer programs were run on the IBM 360 computer using the 
FORTRAN IV language, Since FORTRAN IV does not allow zero subscripts 
we must increase our indices by one and thus take N = 11 to obtain 


ten increments of length ax = Ay = 0.1. Then, T(f,J), (J=1,2,...,11) 


2k+1 


0,J 


corresponds to oy 


corresponds to v » (j=0,1,...,10), and bk Ge tae i) (iene =e eae 
, (1=0,2,...,10); at timem(2ierieat. 
Using this notation, and the approximations (16) and (17), the 


difference equations at the boundaries x = 0 and y = 0 are 


T*(1,J) - Rel) x 2T*(2,J) ~ 21d) . TO nd )o2T) Seth Sth) 
st/2 = ( Ax) ( Ay 


or 


(11b) 2( p41) T*(1,J)-2T*(2,J) = T(1,J-1)+2(P -1)T(1,J)+T(1,J+1) 
for d= 20e i oe ,»N-1) 


an YP +1)T*(1,1)-2T* (2,1) = 2T(1,2)+2¢ -1)T(1,1) 


tome = 1, 
and 
T(T,1)-T#(1,1) | T#(I-1,1)-2T*(1,1)4+T*(I41,1) | 2T(1,2)-27(1,1), 
at/2 re 2 a 2 
( Ox) (Ay) 
Or 


(12b) 2( +1)T(I,1)-2T(1,2) = T*(I-1,1)+2(f -1)T*(1,1)+T*(I+1 ,1) 
hor (l=2,3,. eee) 


2(~p +1)T(1,1) =2T(1,2) = 2T*(2,1) + 2(\e -1) Te Gie) 
for l=l 


Now using equations (lla), (11b), (12a), and (12b) we can write 
out the systems of simultaneous equations which are to be solved, 
Theses are: 


Implicit in the x direction 
Z (Pri) TT) - 2T* (2) = D(l) 
~ TU, T+ Zen TAI) - T*(4,o) 


ee = D(2) 
EC) Ee) a) = D(I) 
~T*W3,T) RCP T(N-2,T)- THIN) ote) 
—T*(N-2, 5) + 20Pt)) TN-1, 7) = D(N-1) 


oD 


where D(I) © T(I,J-1) + 2(-1)T(I,J) + T(I,J+1) 
Cres Ee ) (J=2,3,...,N-1), 


and D(I) = 2T(I,2) + eee Tee 1) 
(I=1,2,...,N-1) 


Implicit in the y direction 


Z(P+0T 1) 27 Cie) = D(1) 

-TC, (+ acre Cr, Om 1,3) = D(2) 

ieee (41) Ce a) T(r, Te) = «D@) 
AT (LN-3)420PH 1 T(L eR N-2)- T(r. ‘" a = = (Nez) 
-T(I,N-2)t2(Pti)T(L,N-I) = D(N-1), 


where D(J) = T*(I-1,J) + 2(P-1)T*(I,J) + T+(I+1,J) 
=e 2 eeg lk) (Em 2,3,0eeeNa1) 
and = D(J) = 2T*(2,5) + 2(P-1)T*(1,d) 
(J = 1,2,.,.,N-1) I=1, 

If we call the coefficients of the main diagonal terms B(I) 
and the coefficients of the lower diagonal and upper diagonal terms 
A(I) and C(I) respectively, then 

A(I) = -1, (I = 2,3,...,N-1) C(I) = -1, (I =1,2,...,N-2) 

B(I) = 2(p +1), (I =1,2,..,,N-l), In addition, let F = 2(y-1). 
The algorithm for solving these systems is then 

A= BU), Sid = DO aan) 

fo(1)= B(1)- AC) CIB cr-1) ¥(1)= [D(z)- AC) 81-1) 82) 
(S25 3 eat-P) 
ONE) = CN) ar A(T) = C= C T(r 1) 0) 


HeUN-1) =—eney) T= rn - CCL)T Citi) (82) 
(I = N-2,N-3,...,1) 
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The flow charts for the solution program ({ADMT) and the algo- 
rithm subroutine (TDIAG) are given in Figures 1 and 2, Twelve time- 
STEPS, BAte=-07 001" 0,00), 0.001. 100002; 07005,50701 770.01, 0,02, 
0.05, 0.05, 0.05, and 0,05, were used consecutively to obtain solv- 

tions at times 0,001, 0,002, 0,003, 0.005, 0,01, 0,02, 0,03, 0,05, 
0 Ol 5 Oana O25. 
The analytic solution for the problem was obtained by the 


separation of variables technique, and is 
elteee samt) Et 

Ulx.y.t) = (16/1? \> > ene Cos (em-EX Cos(em- MJ 
A program was written to obtain the analytic solution at the grid 
points (i ax, — (i,j =1,Z,...,N-l1), for the times 0,005, 0,01, 
0,05, and 0,250, This solution is called ANSOL(I,J), and includes 
all terms in the double summation for n,m = 1 to 20, The values 
of ANSOL(I,J) were read into the program IADMT for comparison with 
the values of the numerical solution, TEND(I,J), The difference, 
W(I,J) = ANSOL(I,J) - TEND(I,J), and the percent difference, 

Z(I,J) = 100 W(1I,J)/ANSOL(1,J), were calculated for all points with 
the exception of those points on the boundaries where the temper-= 
ature is held at zero, 

Figure 3 shows a print-out of the numerical solution and analy- 
tic solution for time 0,25, Figure 4 shows the difference and per- 
cent difference for the same time, 

Numerical solutions were obtained for several values of At 
and Ax, all at time 0,005, Table 1 lists the values of the numer- 
ical solution v(x,y,t) and the percent difference. Z(I,J). at the 


point (0.6,0.6), The analytic solution for this point is U(x,y,t) 
= 0.9999. 


of 


Table 1, 


Bee 


0.0005 
0,001 
0,0025 
0,005 
0,0005 
0,001 
0,0025 
0,005 
0,0005 
0, 001 
0,0025 
0,005 


All solution programs were run on the IEM 360 computer, 


Numerical Solution and Percent Difference at Point 


(0,6,0,6) for Different Values of 4t and Ax, 


b 
Pa 


ee SS See Se SS Se ee 
rPeMrerMro;OCOoOCOoOOCOOoCOo 
iA nn tn Aan MD ODD DO DO 

LA Ln kn ln 


( Ax)* 


0.00062 
0, 00062 
0,00062 
0, 00062 
0,0025 
0,0025 
0,0025 
0, 0025 
0,01 
Ovo 
0,01 
0,01 


At + 


0, 00112 
0,00162 
0,00312 
0,00562 
0, 0030 
0,0035 
0,0050 
0,0075 
0,0105 
0,0110 
0,0125 
0,0150 


AX 


2 


VAX 


0.9996 
0.9995 
Ooo 
0.9981 
0.9994 
0.9994 
0.9990 
0.9982 
0.9975 
Wet 
07973 
0.9965 


i 


Z(I J) 


OF 633 
0,044 
0, 088 
0,180 
0,055 
0,055 
0,099 
0,170 
0, 240 
0, 240 
0, 260 
0, 340 


This 


computer uses seven digits to the right of the decimal point in a 


decimal number, 


off ereor C x 107°, where C is a positive constant, 


Therefore, we can write as a bound for the round- 


(See McCracken 


and Dorm [10] .) C is determined by the number of arithmetic 


operations performed and by how the round-off accumulates during 


these operations, 


large, 


lead us to believe that random cancellation of round-off has occurred, 


thus keeping C small, 


It is therefore possible for C to become quite 


However, the satisfactory results obtained in all solutions 


Furthermore, since the discretization error 


depends on At and ( ax)”, and the smallest of these values was 


©. 2 fia?” we shall neglect the round-off error and examine only 


the discretization error, 


From Douglas 


quantity C, (At) + em tte)™ where C 


1 and C 


Z 
and depend on the least upper bounds of the derivatives U 


[9] , the discretization error is bounded by the 
are positive constants 


and 


U . We can write this bound as K [at + (a x)* | » where K = maxC, »C 


KAKA 


Thus we say that the discretization error is of 0 | ot + (ax) ] ° 


For convenience, we have used At x 10? and ( hx)* x 10° when 
working with the data from Table 1, This data was used to plot 
percent error versus { at + (a x)* | De 10?; the result is shown in 
Figure 5, The line shown in this figure was obtained by a least 
squares fit, Although the data appears to fit a straight line, we 
can not take this as verification that the discretization error is, 
in fact, of 0 [ at + (ax)? | . To obtain verification we would 
need to know, or at least have some estimate of, the derivatives 
Det. and Urye Knowing the analytic solution, it would be possible 
to calculate these derivatives; this was not. done since the objective 
of the paper is to illustrate the use of the IAD method rather than 
to verify predicted results, 

We now compare the work required using the IAD method against 
that required when using the explicit method, We make this com- 
parison for the solution at time 0,25, Although the IAD method is 
stable for any size time-step, we have seen that the discretization 
error depends on At, Therefore, we avoid using very large time- 
steps. For this problem, 0,25 was the final value of t, and solutions 
were obtained for 11 smaller values of t by using values of At rang- 
ing from 0,001 to 0.05, To find the number of arithmetic operations 
required for the IAD method in this problem, we take as our start- 
ing point the solution of equation (lla), which gives the temper- 
atures implicit in the x direction, As we have seen, this equation 
results in Nel = 10 equations in N-l unknowns, The algorithm used 
to solve these equations involves three expressions, each of which 
required three arithmetic operations (including addition and subtrac- 


tion), These expressions were formed a total of N-2 times, so that 


of 


9(N-2) arithmetic operations are required for the solution of each 
system of N-l equations, To obtain the temperature at every grid 
point at the end of a half time-step, equation (lla) must be used N-1 
times, Thus for a half time-step, 9(N-2)+(N-1) operations are requir- 
ed, Now, the solution of equation (12a), which gives the tempera- 
tures implicit in the y direction, requires the same number of oper- 
ations for a half time-step. Thus the solution for a complete time- 
step requires 2x9(N-2)*(N-l) operations, As stated, a solution at 
time 0,25 was obtained as the twelfth complete time-step, so that the 
total number of arithmetic operations required using the IAD method 
was 

12x2x9(N-2)°(N-1) = 19,440 
Use of the explicit equation (9) requires that Ot 2 pieaeehts = 
0.0025, This means that the solution must be stepped through 100 
time-steps to reach time 0,25, The solution of equation (9) in- 
volves 7 arithmetic operations, and this equation must be solved at 
each of the (N-1)° grid points, A total of 100x7(N-1)° = 70,000 
operations are required for a solution at time 0,25, Thus we see 
that the IAD method requires two-thirds less work than the explicit 
method for this problem, This problem involved 100 grid points, As 
the number of grid points increases, it becomes even more advantageous 
to use the IAD method rather than the explicit method, Note also 
that the maximum time-step used for this problem was 0,05, which is 
20 times larger than the maximum allowable time-step for the explicit 


method, 
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uxample 2, Solution of a Steady State Problem Using the IAD Method, 


Problem: To find the steady state temperature distribution at 
the face of an infinite prism with edge-length one. 


The temperature distribution for this infinite prism (See 
Figure 26) is governed by the elliptic partial differential equation 
VU 06y) + U Gay) = 0, 
with boundary conditions 
U(O,y) = U(x,0) = 0, 
Wee) = 1 
U.y) = -hu(l,y). 
The last boundary condition means that the boundary x = 1 is a 
partially-conducting boundary, For convenience we take the value of 
the constant h to be one, 
In solving this problem by the IAD method, we treat it as the 
limiting case of the time-dependent problem Ue Uy = U,. This 


allows us to use the same difference equations, (lla) and (12a), 


that we used in Example 1. 


J 





Ree. og x 


U=O 
figure 26, Temperatures at the Face of an Infinite Prism, 





vie divide the region shown in Figure 26 in exactly the same way 
in which we divided the region for Example 1, Thus the resulting 


grid is identical to the grid shown in Figure 25. 
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THe IAD method, as applied to this steady state problem, can 
be considered as a process of iteration. In such a process we use 
a starting approximation to the solution, and then attempt to i: 
this approximation, We call each attempt at improvement an iteration. 
Thus we regard the use of both equations (lla) and (12a) as one com- 
olete iteration, and we regard the quantity re appearing in these 
equations as an iteration parameter, We now wish to determine the 
optimum values of ~ » i,e., the value or values of e which will give 
us the most improvement in our starting approximation after the small- 
est number of wreratvons. 

Recall from the discussion of the convergence of the IAD method 
that, for a rectangular region, it is possible to predict optimum values 
a cir In that discussion, the effect of one complete iteration was 
represented by i = P( Pi JE, 
where PC py) is a matrix operator and k, is the difference between 
the true solution, U(x,y), and the numerical solution, v (x.y), 
obtained after the kth iteration, We call EY the error at the kth 
iteration, We saw previously that we were assured convergence (for 
an infinite number of iterations) if the eigenvalues of the operator 
P were all less than one in absolute value, These eigenvalues are 


given by, 
x 
(15a) (2s 4 sin PSe (p-4sin’g ot) 
z BX Le 
(p+ 4sin'p $2)(p + 4 Sintq 4) 
(The difference between this expression and (15) is due to the fact 


that the region for this problem has side-length one), 
From (15a) we see that each of the N-l eigenvalues of P can be 


reduced to zero on some one of Nel iterations if we use values of 


ie given by 
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(8) P, = 4 ein pi /2Ne (p = 1,2, «000N-)). 

Thus, theoretically we should be able to reduce the error to zero 
after N-l = 10 iterations, Of course, complete error reduction in 
10 iterations couid happen only when using an exact numerical pro- 
cedure, Due to round-off error we do not have an exact procedure, 
so we can only hope that error will be reduced to a very small, and, 
therefore, acceptable level, Moreover, since in the IAD method, one 
complete iteration involves the solution of both equations (lla) 
and (12a), we require 20 intermediate iterations for maximum error 
reduction, 

Using equation (18), the values of io and the associated values 
of At (from at = (ax)"/p) were calculated and are listed in 


Table 2, 


Table 2. Calculated values of the iteration parameter ss 
Pp. £p At 


1 0, 02462 0, 40612 
2 0, 21799 0, 04587 
3 0, 58579 0,01707 
4 1,09202 0,00916 
5 1, 68713 0.00593 
6 22311287 0,00432 
? 2.90798 0, 00344 
8 3.41421 0,00293 
S, 3.78201 0,00264 
10 3.97538 0,00252 


Since this is a steady state problem, At has no real meaning. We 
use 1t here only to emphasize the similarity in the time-dependent 
and steady state ovroblems when using the IAD method, As in the time- 
dependent problem, where we obtained solutions in order of increas- 
ing at, we also use the values of "a given in Table 2 in order of 


increasing At, i.e., we use the larger values of ‘a first, Each 
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value of ad must be used over a complete iteration so that any 
Saruleuear ‘a is used for two intermediate iterations, 

It is possible to reduce the number of iterations by grouping 
values of a From Table 2 it can be seen that for larger p the 
difference between values of ie is small, Therefore, an average 
value of P may be used for certain groups of Po" The eigenvalues 
of P corresponding to the values of p falling within that group, 
while not zero, are sufficiently small so that effective error 
reduction still occurs, For reasonably small N (as in this case), 
it is possible to calculate all values of 2 and set up the groups 
by inspection, A procedure for setting up groups when N is large 


ls given in (7] . For this problem we group e as in Table 3, 


Table 3. Grouped Values of 2 4 


=e avg At = (dx)</ot 
8,9,10 0.00270 3.7037 
BP or 0,00455 297s 
3,4 0.01310 0.7633 
2 0.04587 0.2179 
1 0.40612 0, 0246 


Thus, by grouping values of et we are able to reduce the nun- 
ber of iterations required by one-half with the expectation that the 
results obtained should be nearly as good as those obtained when 
using all calculated values of a Douglas and Peaceman fi) give 
examples in which simple values of id are used, They start with a 
value between 4 and 10, and then halve this value repeatedly to 
obtain additional values of YP. We shall obtain solutions to this 
problem using all three approaches, 

For a starting approximation, corresponding to the initial 


condition in the unsteady state problem, we use 


dyly 


(19) Ulx,y) = xy/(l-y +x). 
This expression satisfies the boundary conditions at all but the 
conducting boundary, 
As in Example 1, we can use equations (lla) and (12a) at all 
interior grid points and at the boundaries where the values of 
U(x,y) are known, At the conducting boundary, we may use the same 


t 


= Again we assume a fictitious point, a distance Ax outside the 


approximations for U, and “Sa but must find a new approximation for 


boundary, at which the temperature is U(x+ax,y), and as before we 
write the Taylor expansions for U(x+ax,y) and U(x- Ax,y) in terms 
of temperatures U(x,y) on the boundary, where (x,y) represents the 
point [nN bx,jdoy] . Combining these expressions, we have, as before 
(20) U(x- ax,y) + U(x+ ax,y) = 2U(x,y) + (ax) Uy) + NOE. 
However, we no longer have the symmetry which allows us to say that 
U(x= ax,y) = U(xt+ax,y). We must use the boundary condition to 
evaluate U(x+4x,y). The method here is one due to Schmidt as 
given in Carslaw and Jaeger ce . We use a symmetric difference 
for the normal derivative at the boundary to get 
UL(xy) =e ax)| U(x AXy) - U(x- ax,y)| + Wee e 
Solving for U(x+ ax,y), and using the boundary condition 
U(x,y) = -U(x,y) 
we have 
U(x+ ax,y) = Ul(x-ax,y) - 2AxU(x,y) + 0( Boe 
Substituting in (20), we have, finally, 
(21) UL. xy) =e) ( Igoe re QX,y) = (1+ mOUCEy) + O( Ax)*, 
In the computer program IADM we use TST(I,J) to denote the 


approximation to Uix,y) at the end of each intermediate iteration 


a5 


and TEND(I,J) to denote the approximation at the end of a complete 
iteration, In the equations that follow we abbreviate TST(I,J) by 
T*(I,J) and TEND(I,J) by T(I,J), Using this notation and the — 
approximation (21), the difference equations at the conducting 
boundary are: 

Implicit in the x direction 


T*(N,J) - TUN) == 9 2T*(N-1 J) = 201+ Ax)T* (Nd 
At ee 


TIN J-1l) - 2T(N,d) + TIN Jt 


+ 
(ay)* 


or 
(lle) -2T*(N-1,0) + [p(c) + 2@+ax)] T#(N,J) 
= TN,J-1) + [P(k) - 2] 1(N,J) + T(N,J41), 
where 
P(k) = (ax)*/at, = (Ay)*/ ot, 
and k denotes the iteration number, 


implicit in the y direction ; 


TIN,J) = THIN J) = 3 2T*(N-1 J) = 2014+ 4x) TFN 
At OSs F 
n T(N J-l) - 2T(N,J) + T(N,J+1), 
2 
(Ay) 
or 


(12c) -T(N,J-1) + [(k) + 2] T(N,0) - T0N,041) 
= 2me(Nel,J) + [ (kc) - 201+ ax)) P#(N,0). 
We empahsize here that the value a must be used for a complete 
iteration before introducing the value Plkt1). 
Now, using equations (lla), (12a), (11c), and (12c) we write 


out the systems of equations to be solved: 
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Implicit in the x direction 


(2+ yK)] T 2,7) - T°(3,0) = D(2) 

=e "(2, Tyla POT, q)- T* (+, J) 7 oe) SY, 

-T#r4,J) + [2+ pox] T% 2,3) ae r+, d) : ? » D(I) 
~T*(n-2,J) + [2+ PCK)] Tn, 9) at (Nis) = D(N-1) 


—T*(n-1, 7) + [ CK) #21 40) ] (NJ) = D(N) 
(J=2,3,... ,N=1) 
where D(I) = T(I,J-1) + [P(k) - 2| (I,J) + T(1,J41) 
(T=2535... NN) (J=2,5,...,N-1) 
As in Example 1, the coefficients are A(I), B(I) and C(I), with 
B(I) = 2 + P(k), (1=2,3,...,Ne1), BIN) = P(k) + 2(1+ 4x) 
ACI) = -1, (lee es N-1), ACN) = = 
Coe 50 2,oyeet-1) 


Implicit in the y direction 


(2+ (kK)]T (1,2) - T(1,3) = D(2) 

See 2) + [at POT CZ, 3) B)- Tat) — = D(3) 

go J-| Dae PONT ST 5) = T(t, T+) - _ we 
-T(1,p N- 3)+[24POITCIN N-2) — _ Tu, N-1) me Da2) 
-TCI,N-2) + [2+ PCR] TCL, N-1) = D(N-1) 


(1=2 Sate 
where D(J) = T(I-1,J) + [y(k) - 2] T#(1,5) + T*(I4,J) 
(T=2, 3320800) (J=2,3,... »N-2) 
D(NeL) = T+(I-1,N-1) + (yc) - 2] T#(Z,N-L) + T#(141,N+2) 
+ T(I,N) 
(1=2,3,...,Ne1) 
and D(J) = 2P*(Nel,J) + [(k) = 2(1+ax)] T#(N,U) 


I=N | (J=2535 60% ,N=2) 


We 


D(N-1) = 2T*(Ne1,N-1) + [P(k) = 201+ 0x)] T*(N,Nel) + T(N,N). 
we note that D(N-1) contains the term T(I,N) which has been trans- 
ferred from the left-hand side since it is a known quantity, by 
virtue of the boundary condition U(x,l1) = 1. 


The coefficients are written 


B(I) = 2+ p(k), Cl=Z oe ° »N-1) 
Cia, ale (153.4, age) 
CC = -l, Gee AR eC) ° 


We use the algorithm given in Example 1 to solve these two 
systems of equations, The flow charts for the solution program 
IiDM and the algorithm subroutine TDIAG are given in Figures 6 and 2, 


The analytic solution to this problem was obtained from Churchill 


[13] and is 
An An smnhany , 
2a U(x, = 
a2 VRE A>, Kn sinhX, XK, be or 
where AA = ]l =-cosX 


1 + cose es 


and Ky» X53... are the roots of Tan = =X, 

5 program was written to obtain the analytic solution at the grid 
points (i Ax,jay), (i,j=2,3,...,N), and this solution is called 
V(I,J). To obtain V(I,J), the series in (22) was summed until a 
term was reached for which the absolute value of that term, divided 
by the sum of previous terms, was less than 107°, The summation was 
terminated upon reaching such a term, The values of V(I,d) were 
read into the program IADM to allow comparison with the numerical 


solution, TEND(I,J). IA\DM calculates the difference, 


N(I,J) = V(I,J) - TEND(I,J) and the percent difference, 


Z2(I,J) = 100 w(I,J)/V(I,J) for all grid points except at the 
known boundaries, x = y=0andye=l,. 

The first solution was obtained using the five grouped values 
of \p listed in Table 3, Figures 7 through 16 give the results for 
all five iterations. In these figures as well as in thosewhich 
follow, the odd numbered figures show the numerical and analytic 
solutions for a particular iteration while the even numbered figures 
show the difference and percent difference, Results of all five 
iterations are shown, so that the effect of each value of yo and 
the process of error reduction can be seen, We note that error 
reduction takes place diagonally from the larger grid points to the 
smaller, This is because the larger values of ‘i which are used 
first correspond to the larger eigenvalues of the operator P, and 
these eigenvalues in turn correspond to larger values of x and y. 

The second solution was obtained using the ten values of a 
listed in Table 2, Results of the tenth iteration are shown in 
Figures 17 and 18, Comparison of Figures 16 and 18 shows that the 
results of using the grouped values of a in five iterations are 
essentially the same as the results of using all calculated values 
ot 2 in ten iterations. In both cases errors range from approxi- 
nately 0,05 percent to 2.5 percent. We therefore conclude that by 
grouping values of ie we can reduce the number of required iterations 
Significantly, 

The third solution was obtained using the five simple values 
of Pp: 4,2,1,0.5, and 0,25, Results are shown in Figures 19 and 20, 
It can be seen that error reduction was less than in the solution 


obtained by using five grouped values of p. However, these results 


fe 


are considered to be satisfactory when the resulting simplification 
of the problem is taken into account, 

We use the results of the second solution (Figure 16) to fice 
a comparison of the work requirement for the IAD method and the work 
requirement for another iteration method. In the IAD method, the 
equation U(x,y) = xy/(1 - y + x) was used to obtain starting values 
at each grid point, Comparison of these values with the analytic 
solution shows a maximum difference (maximum initial error) of 
approximately 0.10. Comparison of the numerical solution for the 
last iteration and the analytic solution shows a maximum difference 
(maximum final error) of 0.0041, Thus the reduction of error is on 
the order of xi 07“, In Example 1 we found that a complete time- 
step using the I4D method required 2x9(N-2)°(N-1) arithmetic oper- 
ations (including addition and subtraction). For this problem a 
complete time-step corresponds to a complete iteration, and since 
five complete iterations were used, the total number of operations 
required for an error reduction of hx107* was 5x2x9(N=-2)«(N-1) 
= 8100, Peaceman and Rachford [7] give the work requirement for 
the extrapolated Liebmann iteration method of solution to Laplace's 
equation as 5 Y(N-1)? for an error reduction of 107» ; ~This,-tor 
an error reduction of non. the number of arithmetic operations by 
the extrapolated Liebmann method is 5x2(N-1)- = 10,000, Thus the IAD 
method requires fewer operations for a larger error reduction than 
the extrapolated Liebnann method, considered to be one of the more 


efficient methods of iteration, 
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uxample 3, Solution of a Time-Dependent Problem in a Non-rectang- 
ular Region Using the IAD Method, 


Problem: To find the temperature with time for a quarter section 
of the face of an infinite cylinder with radius one. 

The temperature at the face of this infinite cylinder (See 
Figure 27) is governed by the following parabolic partial differ- 


ential equation 


with boundary and initial conditions 
u(r,0,t) = u(r, 7/2,t) =1, u(1,6,t) =1 and u(r,0,0) = 0 


As before, we choose the time unit such that the diffusivity k is one, 


Figure 27, Temperatures at the Face of Infinite linder, 





We divide the region shown in Figure 27 into N increments of 
length ar = 1/N in the r direction, and into N angles of magnitude 


4@ = 11/2N in the © direction. The resulting grid is shown in Figure 


(NOV, No @) 





(/av, NAQ) (Nav, )49) 


=e . 


D 
x (AY, 148) 
Y st V6 
ae 
CLAV,0) (NQV,0) 
Figure 28, Grid for the Region Shown in Figure 27. 
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Let v(r,0,t) be the difference approximation to u(r,@,t), We 


then use the following difference approximations for the partial 


derivatives Uys We? Way? and Ugg: 


VOY, 6, t+at)- VY, Q,) + O (jt) 





a Ae 
Mh = Viv+ av, 6,t)- V(v-Av, 0,¢) 3 Otav)* 
ZO a 
x Viv-ave,t)- Zviv,et + ViVv+av,o¢ 
Vey = uc a a Ovwyv)* 


VV, 8-40, t} - a tV(Vj8400)t) O(0e)* 


Vee 


Now consider the two terms us +(i/r)u,, of (23), They are approxi- 


mated by 
Viv+avie,t)- VVV-aye,t) 4 VW-av,0,t)-Z2ZV(G et) +VVtav, Ot) 


(I~ FV (v-ov,6,t) — 2vW, 6, t) + (14 A¥ 
In Be)¥ ee ON ai 79it) + Ofav)* 





The complete difference equation approximating (23) is 
v(V)0,t+at) - V(V,e,t) 
AT/Z 
-]C- Sv &-o%e, t)-42V,0,t) + C+ 2 )viviave,t)| 
a 
avy se 2VM0,t) + ees) | 
~ 7 t(ae)* ee ~=~=CS—sti‘“Cs—sSN 


= [ O (at) + O(ay)” + O (2e)* | =e) |: 


(24) 





Now consider (24) as an operator L (v) = 0, If L. operates 


on u(r,8@,t), we have 


U(V,O,t ot) - UG O,t) hoy 2 
L (a) = an sce + Ure + OC4t) 
_ [a- %)u@aye,4)-2uM%9,t) + 41 )Jucraye, t)] 
(ery 


= Be 080)" pery + Olav}s_ UMe-49 t]- 2d 2, t)+ 4 Yotsg,t) 
4 v*(Ae)* 
~ CN Geen. + Oe) = O 


Te a 





De 


where the tilde indicates that the derivatives are to be taken at 
certain points in the region 0<r<1, 0<6 <7/2, O<t ST. 
We call z(r,0,t) = v(r,0,t) - u(r,0,t) the discretization error, 
Then since L. is a linear operator, we have 
Lv) -Li(u) = Liv - x), 
or 


~ ~ tw ca of q 
L (2) = Sele 2 1 EY (av) Lay ts = Ugesy t+ Oat) + O@v)"t OG) 


Call the right hand side of this equation w(r,®,t), Thus we see 
that z(r,6,t) is a solution to the difference equation 
L (2) = w(r,8,t) 

Since (24) represents a convergent difference scheme, we have 
that 

lusb [zl = lub iwl, 
or 

(25) meee ze C,4t + C.( ar)* + C( ae)”, 

nerr 


where C,, C, and ©, depend on the least upper bounds of Uzppe U 


and respectively, In addition, C, and C 


“9800 2 3 
fore, we Should expect a discretization error of 0 [ at + (ar)*+ (ney'|. 


depend on r, There-= 


We now digress somewhat to consider the IAD method as applied 


to the related elliptic PDE (u, +(1/r)u,, +(1/r*) up, = 0) for this 


r 


non-rectangular region. If U(r,@) represents a vector solution to 


this equation, then the difference analog to (24) is 


~ on Lab +A°] u(r,e) = 0 


D9 


where Ne and A’ are square matrices such that 


aPu(r,e) = - (i= SY) u(r-ar,@) + 2U(r,8) - (1 + ¥)U(r+ ar,e) 
ay) 
A’U(r,0) = - 2O, [ver.e- 46) = 2U(r,8) + U(r ,6+26) | ; 


For the simple case of N = 4 (a net with 9 interior points) we 


find that 
2 -3/2 0 0 0 0 0 0 0 
-3/4 2 -5/4 0 o oOo 0 0 0 
0 -5/6 2 0 0 0 0 0 0 
C 0 0 -’)  ® 0 0 
a 0 0 0 -3/4 2 -5/4 0O 0 0 
ee 0 0 oO oO 5/6 2 0 +O 0 
oO oOo oO Oo oO 0 2 -3/2 0 
oO oOo 0 0 860 0 -3/4 2 =5/4 
0 0 0 0 0 0 0 -5/6 2 
and AY is similar to a matrix B such that 
2 el 0 0 0 0 0 0 0 
-1/4 1/2 -1/4 0 0 0 0 0 0 
0 -1/9 2/9 0 0 0 0 0 0 
0 0 0 2 -l Oo Oo 0 0 
B = 0 0 oO -1/4 1/2 -1/4 0 0 0 
0 oO oO oO -1/9 2/9 0 OO 0 
0 0 0 0 0 0 2 -l 0 
0 0 0 0 0 QO -1/4 1/2 -1/4 
0 0 0 0 0 0 0 -1/9 2/9 
V 


For this case, although fi and A are non-singular, neither is 
symmetric, and a calculation shows that ADyv ¢ sy Thus, for 
this non-rectangular region, ae and A’ do not share a common eligen- 
vector basis, and it would therefore not be possible by the method 
used in Example 2 to predict values of the iteration parameter for 
optimum convergence rate, It would be necessary to use simple 
values for this iteration parameter (See Example 2). This is con- 
sidered to be the most serious disadvantage of the IAD method, 
However, for the parabolic problem under consideration, we shall 
see that entirely satisfactory results are possible even though 


the region is non-rectangular, 
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We now return to equation (24), Call ve j the approximation 
9 


to u(r,0,t) at the grid point (iar, j40, kat). Let P= (ar)*/at 
and X = 4@; also note that r=iar. Then equation (24) yields the 
"working equations" analogous to equations (lla) and (12a), These 
are: 


Implicit aagthe, xr direction 








(26) = Cea oh f LCP) Ey zz) Venn 
= ax Vey + 2 (eee) Vey + re Mest 
Implicit in the @ direction 
(27) = a ot 2(+. oa PVE — ae _ shina 
= dan + BPO Cora Men 


Since the coefficients in these equations depend on i, this 
will add some complexity to the computer program, However, since 
the temperature is given at all boundaries, we need find no special 
approximations at any boundary, as was necessary in Examples 1 and 2, 
We use the same notation used in the previous two examples when 
writing out the two systems of simultaneous equations which result 
from Kquations (26) and (27), These are : 


Implicit in the r direction 


i 

o 
oo 

Nh 
Ww 


2(p+1)T#(2,J) - (141/2)T*(3,d) 


(1-1/4) T*(2,5)+2(9 +1 )T*(3,d)-(141/4)T*(4, J) = D(3) 


-1(1-1/2(N-3) )T* (N=3 ,J)+2QP+1 )T* (N-2,J)-(141/2(N-3) )T*(N-1,J)=D(N-2) 


~(1-1/2(N-2) )T* (N-2,J)4+2(o9 +1) T*(N-1, J) = D(N-1) 
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where 


D(I) = (FE ; Eye = = Ge ie = ee ae .) T(2,J) * yea aie Gr J+1) 
WIEZ 3 yee Nae (T=2 ,3,0se9hl=L) 
and 
D(2) = D(2) + 0,5; D(N-L) = D(N-1) + 1 + 1/@(N-2)) 


The coefficients needed for the algorithm subroutine, TDIAG, are 
B(I) = 2( P+) (122335 eee NL) 
(= =e), 20 ee ae 
C(I) = -(141/2(I-1)) (1=2,3,...,Ne2). 


Implicit in the 8 direction 


2|¥ nonkae ENE Cl S) aE NED) 
Tenge t+ ALP are mel a < ral pr = D(3) 
cere ET) + [Ft ere (TE I)- Cee TD It) = D(J) 
aapege TUN) 42 [Y + pe] TC Ne) TCE Nt) = D(we2) 


sage TG N- +] pt ote TER? T(ZN-1) = D(N-1), 


where D(J) = terion a *(I-1,J)+ 2(P-1)T*(I, 7) + [1+ 


(T=2508... ,N=1) ) (J=2,3,.--,Ne1), 


Tray sh T HHT) 


and D(2) = D(2) + Se D(w-1) = D(N-1) + TTX 


with coefficients 
B(J) = 2¢P41/(1-1)*x*)) (d=2,3,... Ne), (1e2,3,... MeL) 


- | 
A(J) = (I-1)* «* (J=3,4,...,N=1) ) (1=2,3.... »N-1) 


C(J) = A(J) (J=2,3,...,N=2) C7 ee 


56 


Note that in these equations, each coefficient has I reduced 
by one, This is a consequence of having to avoid zero subscripts, 
The flow charts for the solution program, IADMC, and the 
algorithm subroutine, TDIAG, are given in Figures 21 and 2, 

The analytic solution is given in Jaeger (14] for the points 
r = 0,25, 0.5 and 0.75; @ = 15°, 30° and 45°(the temperature is 
symmetric about the 45 degree line) at times between 0,01 and 0.2, 
These temperatures are given to three significant figures, but the 
author states that calculations were made to four significant 
figures so that the third figure is considered to be correct. 

Figure 22 shows a printout of the intermediate numerical sol- 
ution, TST(I,J), and the numerical solution, TEND(I,J), for time 
0.1. For this solution we used ar =1/12, 40 = ™/24 and at = 0,005, 


Table 4 is a comparison of the analytic and numerical solutions at 


tien. 
Table 4. Analytic and Numerical Solutions for Time 0,1 
r—»> On2> O25 OF ee) 


anal. Sole ilumeeool. | Anal, Sol. |Niums Sol, PAwal, Sol.| Num. Soll, 


{, 


15°} 0.979 0.9788 0,963 0.9631 0.957 0.9574 
30°| 0,946 0.9453 0.906 0.9053 0.891 0.8907 
45°} 0,952 0.9518 0.917 0.9166 0.904 0.9037 


Nine different solutions were obtained using combinations of 
ar, 40 = 1/12, 1/24; 1/24, 1/48; 1/36, 1/72 and at = 0.005, 
0.0025, 0,00125 for time 0,02. Table 5 shows these results at point 
(0.5, ES). The difference listed in this table is the difference 
between the numerical solution and the analytic solution at this 


point, The analytic solution is 0.168, 


oy 


Table 5, Results +f solutions using various 4r, 486, and St 
at _ the point (0,5, 45°) 
At + Num, 

At Or IN: Ar)*+( 06 Soln Difference 
0,00125 0,0277 0, 043 0, 00393 0,167 -0,001 
0,0025 0.0277 0.0436 0,00517 0,168 0 
0,005 0,0277 0,0436 0, 00767 0,171 0,003 
0.00125 0,0416 0, 0654 0.00728 0,168 0 
0,0025 0, 0416 0,0654 0,00853 0,169 0,001 
0,005 0,0416 0.0654 0, 01103 0,172 0,004 
0,00125 0, 0833 0,1309 0,01345 0,174 0,006 
070025 0, 0833 0,1309 0,01470 0,175 0,007 
0,005 OSs 0,1309 0,01720 On76 0,008 


Figure 23 is a plot of the absolute difference versus At + 


Ceoe + (a6)* using the data given in Table 5, The line shown in 
this figure was obtained by least squares and has the equation 
0,626x - 2,88 (we have scaled all values by a factor of 10°, for 
convenience), Here, as in Example 1, the data appears to fit the 
straight line, but, as previously stated, we cannot take this as 
verification that the discretization error is of 0 [ Ot + Cae 
+ (ae)*] since we do not have a predicted value for K, The 
analytic solution as given in faa] is 


0° —_ 
u(r,0,t) = I-22 Sean thin J Sg Msn) e Jig (Memt dt 


SBA me (F/G sm))” 
where J. is the Bessel Function of order s, (s = 2(2n +1)8, 
n=0,1,...), and re: (m= 1,2,...) are the positive roots of 
J_(K) = 0, Thus we see that it would be very difficult to obtain 
the derivatives User Une? and Upgeee needed to predict a value for K, 

The complexity of the above equation also leads us to conclude 
that the IAD method yields a simpler program than one involving the 


analytic solution, 


Another advantage of the IAD method in this example is that, 
by solving the equation in polar coordinates, we avoid the necessity 
of interpolation at the circular boundary. Such interpolation be- 
comes necessary when placing a rectangular net on a circular region, 


Since not all of the grid points will fall on the curved boundary, 


D7 


CONCLUSIONS 

From the numerical results of the three examples, we have 
seen that accurate results are possible when using the IAD method, 
The results for the parabolic problem in the circular region suggest 
the applicability of the method to many types of problems, including 
fluid flow and electric potential. When the differential equation 
is written in polar coordinates,a possible difficulty arises when 
dealing with derivative-type boundary conditions at the origin, 
Since the differential equation has a singularity at that point. 
(Temperatures on the boundary were specified in Example 3, thus 
avoiding this problem,) Such a boundary condition would require that 
additional approximations be made, However, the method seems 
ideally suited to problems involving annular regions, As pointed 
out previously, one advantage in using the IAD method on equations 
in volar form is that there is no need for interpolation at curved 
boundaries, A second advantage is that the computer program for 
the IAD methodis probably simpler to write than a program to obtain 
temperatures at the grid points using the analytic solution (if know), 

Since Example 2 involved an elliptic problem over a rectangular 
region, we were able to calculate optimal values for the iteration 
parameter ‘e . For non-rectangular regions the same procedure is 
not applicable; weweu, therefore, be in doubt as to the number of 
iterations required, However, it is vossible, when using simple 
values of p , to gauge the progress by calculating the sum of the 
squares of the errors after each complete iteration to see what the 
effect of a varticular value of id might be, Details for this pro- 


cess as well as other examples of solutions in non-rectangular 
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regions may be found in Douglas and Peaceman (11] » I addition, 
Wachpress [16] suggests a combination of the IAD method and an 
iteration method due to Lanczos for elliptic problems in non-rect- 
angular regions, Douglas [25] ; [27] also gives details for 
extension of the IAD method to mildly nonlinear problems and to 
problems in three space variables, 

Some doubt remains whether the data obtained in Examples l 
and 3 actually satisfies the predicted error bounds, As we have 
seen, the data obtained from the numerical solutions is not sufficient 
in itself to verify predicted results, The whole question of dis- 
cretization error is found to be much more complex thar. was supposed, 
and, in fact, could be the subject of a separate paper. 

The three examples given show that the IAD method involves 
essentially the same equations for both elliptic and parabolic 
problems, The algorithm subroutine used to solve the simultaneous 
systems generated by the IAD equations is exactly the same for both 
problems and requires only a change of input parameters, This fact, 
and the satisfactory results obtained, coupled with the ease of 
application, leads us to conclude that the IAD method has a wide 


range of application in the field of partial differential equations. 


onl 


LO, 


8 oe 
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DIMENSION A,86,C,D 
TST, TEND, TEMP, OT, 
RHO, ANSOL, W, 2 








006 J=I,NA 6 READ 
TEND(I,7) VALUE S 
DOG I= 1,NA = 1,0 OF OT 






TAU = 
TAUs DOT Ct) 





XINC @XINC 
OT CT) 










B(i)=2,0*RHO+1.C -1,0 
B(J)= BC) 


CCr}= - 1.0 










F= 2.0*RHO-10 


O(1)= 2.0 * TEND (I,2) 
+ F * TEND (I, !) 


D(I)= TEND(I,J=1) + 
FATEND(L,S)+ TEND(L, J+!) | 









p(3)= 2.08 T5T(2,7) 












DO60 JtI,NA 





O(3)= TST(I-i, 5) + 
Fe TST(1,5) + TST(I+I,5) 


FIGURE 1 FLOWCHART FCR PROGRAM IADMT 


OD 





160) 
CONTINUE 


CALL TDIAG 













(NA, A,8,C,0, TEMP) 

















READ VALUES 





TEND(I,J) 
= TEMP(CS) 









OF ANSOL DO!HI0 J=1,NA 














W(L J) =ANSOL (1,5) 











~ TEND (I,J) PRINT PRINT NUMERICAL 
é 
Z(I,J) =100*W(Z,9) TIME, TAU SOLUTION, TEND 
ANSOL(I,J) 







PRINT ANALYTIC 
SOLUTION, ANSOL 


PRINT PERCENT 
DIFFERENCE, & 


OT= 4t 


PRINT 
DIFFERENCE , W 


XINC = ax 





RHO & a 


FIGURE 1 (CCNTINUED) 
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SUBROUTINE TDIAG 
(NC,A,8,C,0, TEMP) 





DIMENSION A,B,C,D, 
] BETA, GAMMA, TEMP 








BETA(I) = BCI) 
- A(I)*C(I-1) 
BE TACI-!) 






BE TA(2) = B(2) 
GAMMA (2) = 
D(2) /BETA(2) 



















GAMMA(I) =[D(CZ) 
- A(I)# GAMMA(1-/)] 
/BETA(L) 


TEMP(NC) 
= GAMMA (NC) 
ND=NC-2 











DO 10 K=1,ND 


TEMPCT) = GAMMA(T) 
— C(I) * TEMPC(I +1) 
BE TACT) 







FIG'RE 2 FLOVCHART FOR ALGCRITHM SUBRCUTINE 
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PERCENT DISCRETIZATION ERROR 


004 








Line ig Y= O22x . 


- 0 ee ae oe ee 


i 
| i 


SS I OS A | A SS Oo eS 
1 


.08 


.06 


FOR TIME 04005 AT PCINT (.6,.6) 


~ “FIGURE Fernie clos nat a + (0x) 
— IAD SOLUTI ! 
| | | 


; 8 /o i2 (4 /é 18 
[at + (ox)*] x 10° 


OF 7 Uxytlyy Uy 
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DO2 J=z\,N 









DIMENSION A,8,C,0 
TEND, TST, TEMP, RHO, 
V,W,z 











00 2 T=i,N 





00 8 J=i,N 









TENOCI,S)=0 DO5 Is1,N 




















Ve Chinas 2 Pra trd x) PRINT TEND(I, J 
X= (XI=1.0)/10.0 TENOC(L,J)= TST (1,7) I=',N Jz 1,N 













8(1)= B41) 
C(t)s-5.0 


6(1) = RHO+ 2,0 
Fs RHO-2,0 

























O(I)< TENDS, Jes) 
+FeTEND (I,J) 
+ TEND(I, J+!) 


A(N)s -2.9 


B(N) =RHO42,0(1,0 + 
KXINC) 


D020 J52,NA 








DO 30 122,N 
NAsN-l 










CALL TOIAG 






TST C1I,T) 
= TEMP(I) 






(N,A,8,C,0, TEMP) 





O(J)= TST (t-1,3) 
+F4TST(I,T) & TST(L+I,T) 


FIGURE 6 FLOWCHART FOR PROGRAM IADM 


fal 





CALL TDIAG 





50) 
CONTINUE 













D(NA) = 
D(NA) +1,0 






(NA,A,B,C,0 TEMP) 














DO 205 J=2,N 
00 205 [=2,N 


PRINT NUMERICAL 
|SOLUTION, TEND 


PRINT PERCENT 
DIFFERENCE, Z 






TENDC(L,S) = 
TEMP) 


D040 J=2,NA 












WZ, J)=V(L,5)-TEND(S, J) 
Z(1,J)=100,.0"w(r,7) 












PRINT RHOCK), 


ITERATION 
NUMBER, K 














PRINT ANALYTIC 
SOLUTION , V 


PRINT DIFFERENCE 
W 





XINC = OX 
RHO = ? 





FIGURE 6 (CONTINUED) 
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PT= 3.14159 
RINC = 1.0/RN 
ALFA = PI/24.0 


DIMENSION A,6,C,D 
TST, TEND, TEMP 


























AS = ALFA¥ALFA 
DT=0.005 
TAUZ0.0 


DOS J=zi,N 
005 [=1,N 


RHO = 
RING € RINC 
OT 







5S} TSTULS)£O TST(II= 1.0 


TENO(I,7)< 0 









TENOC(T, I) E10 







6) TST(LNI=1.0 TSTCUTIZ 1.0 









N 










TEND (I,N)= 1.0 TEND, TZ) = 1.0 












8) TST(IN TI 1.0 
TEND(N,S)= 1.0 





00 10 Jz2,NA 














Afc J-! 
J8(I)= 2.90 
# (RHO + 1.0) 






DO ZO J=2,NA 










=(|.0+1.0/C2,0 AT )) 
AC 


po gO 1=2,NA 
~ (1,0=1.0/(2.0% AT)) 









D(z) = (1.0/CAI#AS))*% TEND (1, J-1) 
+ 2.0 *(RHO=-1,0/(AIMAS))* TEND CI, J) 


+ (1.0/(AZMAS))%TENOCL, T+!) 





(I-1)™ (I-11) 





FIGURE 21 FLOWCHART FOR PROGRAM IADMC 
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30) 


CONTINUE 















0(2) = 042) +0.8 
O(NAS S DINA)+ 1,0 
+ 1.0/(2.04 (RN -1,.0)) 


0040 172,NA 


50} C(k)= -!.0/(AI*AS) 
B(K)=22.0%(1,0/(AT* AS) + RHO) 


AK) = CCK) 


CALL TOIAG 
CNA, A,B,C, 0, TEMP) 















TST(I,J) 
= TEMP(T) 












00 60 Jz2,NA 











|GO 
CONTINUE 





DCI) = (1.0-1.0/(2,.08AI1))* TST (r-i,3) 
+ 2,08(RHO-1,0)* TST (2,7) 


+ C1.041.0/(2,0¥AT1))# TST CI+1,J) 












0(2) 20(2)+ 1.0/(AT¥AS) CALL TDIAG 









90 40 J=2,NA 





(NA,A,6,C,0, TEMP) 





D(NA)= O(NA)+1,0/( AT# AS) 












PRINT 
INTERMEDIATE 
SOLUTION, TST 







TEND(I, J) 
= TEMPCS) 







PRINT 
NUMERICAL 
SOLUTION, TEND 






RING = aY 
ALFA=s 40 
OT= at 
RHO = ~ 


FIGURE 21 (CONTINUED) 
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